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Abstract 

A multivariate dynamic Markov model is formulated to 
describe the possible effect of a generic chemical toxin on a generic 
cell population in an organ. Asymptotic methods (large cell 
population) are used to show that numbers and toxin may be 
jointly Gaussian/normally distributed; the joint stochastic process 
is approximately Ornstein-Uhlenbeck. Application to dose- 
response, and hence risk analysis, is briefly discussed. 



EXECUTIVE SUMMARY 



This paper presents a generic model for the influence upon the cells in an 
c: n caused by dosage by a chemical toxin. The model is presently non-specific 

to either organ or toxin, but provides a starting point for specialization. The 
model is probabilistic/ stochastic without requiring computer simulation (it relies 
on asymptotics appropriate for organs containing a large number of interactive 
cells). 

The organ disfunctionality associated with toxin input that is modeled is the 
(probability distribution of the) number of functionally active cells present: that 
number tends to be reduced by toxin input. Inter-cell signalling is modeled by a 
pairing device; this functionality is due for later elaboration. The effects of cell 
death by necrosis and apoptosis are modeled. The present model omits explicit 
consideration of such issues as spatiality of cells within organs, enzyme function 
and other important features. 

Application of the present model is made to consideration of the shape of a 
dose-response function at low doses. 



AN EXPLORATORY STOCHASTIC MODEL 
FOR TOXIC EFFECTS ON CELLS 

R. L. Carpenter 
D. P. Gaver 
P. A. Jacobs 



1. Introduction 

Application of laboratory toxicology data to environmental and human 
problems of risk assessment almost always requires extrapolation of the data 
from the experimentally used dose regimen to the exposure conditions of 
practical concern, and from the animal species tested to the species of concern 
(usually man). This extrapolation and estimation process is known as chemical 
risk assessment. The risk assessment process has undergone considerable 
evolution in the last ten years, moving from a qualitative basis for decision 
making to an increasingly quantitative basis, and from the use of default 
assumptions to the application of mathematical models as tools upon which to 
base decisions. In the context of determining safe human exposure limits to 
potentially toxic chemicals, there are two sub-tasks to be accomplished: 
estimation of low dose risk in animals and the subsequent conversion of animal 
risk estimates to human risk estimates. 

Dose extrapolation, the first sub-task, can be accomplished either by 
assuming biological response varies linearly with dose, or by using 
physiologically-based pharmacokinetic (PBPK) models. Multi-compartment 
physiological models are formulated using actual tissue volumes from the 
experimental and target species and actual perfusion rates to provide for 
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chemical transport between the compartments. Thus the pharmacokinetics of 
high to low-dose extrapolation become amenable to calculation, and external 
measures of dose can be translated to concentrations in the target organ (internal 
dose). Target-organ chemical concentration may be translated into estimates of 
risk if a mechanistically-based model can be used to relate chemical presence in 
the organ to harmful outcome. Such models exist for carcinogens, in the form of 
the widely-used linearized multistage model. This model is based on the 
generalized concept that chemical alteration of the cellular genes may give rise to 
permanent, heritable changes in the genetic information stored in the cell nucleus 
and lead to phenotypic changes in the altered cells that ultimately cause the 
formation of malignant tumors. 

Unfortunately, no similar models exist for chemical toxic effects other than 
carcinogenesis. The mathematical expression of even the relatively simple 
concept of genetic alteration leading to cancer involves significant simplification 
of biological reality and significant mathematical complexity. This state of affairs 
is exacerbated when one attempts to accurately describe the interactions leading 
to loss of cell function and cell death in tissues of a whole animal. The multi- 
layered control and response systems present in an intact living animal are 
poorly understood and have not been modeled as yet. Nonetheless, these control 
mechanisms defend against the majority of toxic effects observed as a result of 
chemical exposure, and their failure represents most of what is observed as the 
expression of toxicity. 

This paper describes initial formulations of models describing cellular 
response to toxic insult. Our model formulation takes into account aspects of the 
current state of knowledge concerning the control and regulation of cellular 
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properties in tissue. The design of this model is not to provide a detailed 
description of the response of specific organs to toxic insult, but rather to begin 
the process of mathematically describing the recent significant advances in 
knowledge that have been made in cell biology, hormonal regulation, and 
hormetic control mechanisms. Considerable extension of model formulation will 
be needed to apply it to specific organs, since the complex geometry of organ 
architecture is not present. To have included complex, three-dimensional 
relationships between different cell types would have increased the mathematical 
complexity considerably. The current formulation of the model is suitable for 
experimental evaluation in cell culture if the experimental conditions are 
properly chosen. 

2 . Biological Background 

The tissue making up organs almost always involves a complex geometrical 
juxtaposition of several cell types having specific and often overlapping 
functions. This tissue architecture is maintained by the presence of a nonliving 
matrix of proteins collectively known as a basement membrane. Interactions 
between this membrane and the cell are known to be critical to its stability as a 
mature, functioning entity. The whole of the tissue is permeated by branching 
blood vessels, each generation of which is successively smaller; these serve to 
provide a constant, stable milieu in which the cells exist. Nutrients, control 
signals in the form of specific biomolecules, and xenobiotic chemicals are brought 
to the immediate surroundings of the cell by the vasculature in the tissue and 
cellular metabolic products; the products of energy metabolism and cellularly 
derived control biomolecules are removed from the cellular microenvironment 
by the same means. 
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The state of the cell at any time is a reflection of its age, the summation of the 
control chemicals reaching and leaving the cell, the effects of xenobiotic 
chemicals present, if any, and the state of the cells surrounding it. Cellular 
contact with the surrounding cells and basement membrane act to supply 
chemical signals that modulate its activity. A given cell may be (i) nearly 
quiescent, (ii) active biochemically, i.e., producing metabolites of absorbed 
materials for its own use or for export, (iii) in a state of stress due to shortage or 
excess of biochemical molecules, (iv) in a process of programmed cell death 
(apoptosis), (v) in the process of dying from chemical insult (necrosis), or (vi) 
dividing to form progenitor cells in response to a need to replace cells already 
lost; these conditions need not be mutually exclusive at any point in time. Some 
specific cells may alter or completely change their observable characteristics in 
response to chemical signals. The most well known of these cells are the 
pluripotent stem cells of the hemopoetic system. 

3. Prototypic Mathematical Models for Cellular Response to a Chemical 

Toxin 

In this section we propose several simple models to describe the interaction 
between cells and a chemical toxin with which they are in contact. The aim is to 
describe anti-toxin functionality in a parsimonious manner, so certain details are 
omitted that can feature in later models. 

Consider a collection of interacting cells in proximity in an organ such as the 
liver. In particular, these may be the cells that line the sinusoids in the liver, the 
function of which is to remove waste from the blood flowing through. 

Focus first on the cell cycle. For present purposes a given cell may be in one of 
two states: functionally active, i.e. or undergoing mitosis, i.e., mitotic, other 
detailed aspects of the cycle are temporarily ignored. During the functionally 
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active state occupancy time a cell performs its intended function, i.e., that of 
removal of nutrients and wastes from contiguous blood flow. Such a cell has a 
certain natural life time, the duration of which is influenced both by the 
programmed cell death phenomenon apoptosis, but also by the occurrence of an 
insult, possibly directly from an environmental agent such as a toxic chemical, 
but also from a metabolite; the latter cause of death is necrosis. When a cell dies a 
chemical message is sent to a neighboring cell via the intercellular media 
commanding that neighbor to divide, or perform mitosis, i.e., go into S and Cj- 
During the period of a cell's mitosis the original collection of active cells is 
effectively deprived of two members: the dead one, and the neighbor undergoing 
mitosis. Note that a dying cell need not always signal the same neighbor to 
become mitotic, but restriction to that situation is not exceptionally special and 
does lead to quite a simple initial model. 

Model D.l . A Deterministic Model Involving Monogamous Cell Pairs. 

Consider a monogamous cell pair that alternatively is in mitosis and is 
functionally active. Suppose {Xf, i = 1,2, ... } are independently and identically 
distributed positive random variables that represent the time spent in mitosis for 
the pair, and {L]i, L 2 i , i = 1,2, ... } be the lifetimes of the two paired cells when 
they are in the functionally active state. We assume the pair's lifetimes are 
probabilistic, being mutually independently and identically distributed, and are 
re-sampled after each mitosis event period terminates. The pair dies, in effect, 
after lifetime L,- = minfL],-, L 2 O, after which mitosis begins, eventually terminates, 
and a new pair lifetime begins; so the process continues. If Id(D is the indicator 
function that specifies the state of the pair at t, so 

[1 if both cells of pair are in totally differentiated state at time/; 



7 d(0 = 



0 otherwise 



(3.1) 
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then 



p D (t) = P(l D (t) = 1 I Both pair members have just entered the totally 

differentiated state at t = 0} (3.2) 

satisfies a renewal integral equation, see Feller (1966): 

Pd( 0 = 1 - *7.(0 + Jo *b(' - x)G{dx), (3.3) 

where 

h(>) = p { L li 5 ‘Mi S <} = ?{£• = ma x[L v ,L 2i ) £ (} (3.4) 



and 

G(x) = F l (x) *Fx(x) = JqF l (x - y)Fx(dx). (3.5) 

The argument for (3.3): starting with both pair members freshly back from 
mitosis the pair is either, (i), still in the functionally active state at time t, an event 
of probability 1 -Fi(t) --the first right-side term of (3.3) -- or, (ii), the pair has 
cycled, first through its functionally active state and then through mitosis, with 
the pair starting a new life at time x; this is represented by the integral term of 
(3.3). 



Asymptotics 

Renewal theory results, cf. Feller (1966) or Asmussen (1987), tell us that, as 
time gets long, a limit is approached, expressed as 

lfanq,(/).p D . (3.6) 



where E[L] is the mean time of pair residence in the totally differentiated state, 
and E[X] is the mean time of cell pair's absence for mitosis. Thus the productive 
fraction of the time, po, depends only upon the means of the state residence 
times and not upon details of their distributions. Of course different cell types, 
and those with different locations in an organ, might well have different means, 
and hence different pp-values. Moreover, both means could be expected to 
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respond to a toxic chemical, so we would write E[L; T] and E[X; T ] where here T 
stands for level of toxin, past and present. At the moment no attempt is made to 
represent the dependence of the above means on toxin level in a mechanistic 
manner, but that is an ultimate objective. 

Groups of Cells; (Sub) Organs 

Suppose there are C pairs of cells in an organ, with pairs behaving as 
discussed, and independently cycling between states. Then standard limit theory 
suggests that if D(t) represents the number of pairs totally differentiated cells 
present in the organ at time t, 



is approximately normal/Gaussian distributed for large C. Divide numerator 
and denominator of the r h s by C to conclude that as C becomes large the 
concentration of functional cells D(f)/C approaches a constant (in this case) i.e. 



Model D.2 . A Differential Equation Representation of Cell-Toxin Interaction. 

Let D(t) represent the random number of cell pairs that are in total 
differentiated state at time t, and M(t) be the number of (potential) pairs 
undergoing mitosis; D(t) + M(t) = C, a constant representing the number of cell 
pairs in the organ. Suppose the process {D(t), t > 0} is Markov with generator 
defined by 



z D (i) - d(,) -^ D ) 




(3.7) 



PD- 



p{D(t + dt ) = D(t) + l|D(f )} = p{T)(C - D{t))dt + o(dt) 
P{D(f + dt ) = D(f ) - l|D(l )} = A ( T)D(t)dt + o(dt) 



(3.8) 
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Then if T is a constant level of toxicity in the system one can take expectations 
throughout to obtain a linear differential equation for D(t) = E[D(t)] the mean of 

D(t ): 

% = M (T)-C-(rfT) + Z(T))D(t); (3.9) 

at 

the solution of which is 



D(t) = * C f°2 J l - £-(" (7 ‘ )+;l(T)) ' i (3.10) 

w w /t(T) + A(T)V ) 

Hence the long-run expected number of totally differentiated pairs is 



lim D(f) = 

t —>° o 



Cjij T) 

u(T)+m 



(3.11) 



which bears close resemblance to Model 1 in the long run, expression (3.6), with 
= E[X] and A(T,H = E[LJ. For the above model it is again easily seen that (a 
scaled version of) D(t ) is approximately Normal for large C. 

Cell-Toxin Interaction 

The above motivates the following deterministic differential equation 
formulation. We allow (3.9) to become 



o r(t) v p(on o 



(3.12) 



dt w ’ 1+ KT(t) 

The differential equation for T(f), roughly interpretable as the mean 



concentration of toxin present at t, expresses the fact that the increase in toxin 
present at t equals the rate of input of toxin minus the rate of output; the latter 
rate is modeled as a saturated (Michaelis-Menten type) service rate, 
t/T(f)/(l+ x"T(f)), multiplied by the concentration of totally differentiated cells 



in the fluid. In a later section we study a complete probabilistic version of (3.12). 
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It is perhaps natural, but is conjectural, to attribute functional dependence of 
the cell-change rates in a way that tends to reduce the rate of increase of 
functional cells, the “toxin fighters", with toxin concentration. For example 
(only), let 

/r(T) = /v“ |T ; A(T) = A 0 ^ t (3.13) 

where £ and 77 are positive. This choice tends to kill more cells as T increases but 
to increase the number undergoing mitosis since the duration of the mitotic 
period is increased. Another possibility is that the condition of cells emerging 
from mitosis is changed in the presence of a toxin: some such cells are deficient, 
perhaps leading to birth defects, or others become initiated for cancer. We do not 
now model these many options. Supposing r(f) = t, a constant, the long-run toxin 
concentration is the solution of 



- Qi(T) 

D %(T) + A(r) 



(3.14) 



and 



or 



dt _ t £h(t)t 

T ” ^ 1 + kT ~ (l + xf )[/i(f ) + A (f )] 



vCn 0 T 

(l+K'fj^„ + V (i4 ' ,)7 ’] 



(3.15) 



A graph shows that a finite steady-state long-run solution for toxin level, namely 
T, wall exist under certain conditions. Note that in Fig. 1 two solutions to (3.15) 
can exist. 



If 



t < max 
T 



vCfi P T 



(1 + vT) + 
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T 



Figure 1 

It seems clear that the smaller of the two possible solutions is biologically 

plausible, since it is evidently an increasing function of f, whereas the other 

solution decreases. If the above inequality is not satisfied the organ is soon 

overcome and toxin concentration increases indefinitely. 

Model D.3 . A Cell-Age-Specific Differential Equation Model of Cell-Toxin 
Interaction; The Effects of Apoptosis and Necrosis. 

In order to model the effects of programmed cell death or apoptosis, and death 

from insult, called necrosis, it is necessary to represent the cell-aging 

phenomenon. We choose to do this in the context of Markov chains or rate 

equations, simply extending (3.12), by the consideration of two totally 

differentiated cell types: new or young, namely those that have comparatively 

recently been "born", i.e., emerged from mitosis, and old, those that have 

transitioned from the new state ("reached maturity") and are now eligible to die. 

The following rate equations reflect the effects of a toxin in the above process. 

Let D n (t) be the number of new or young cell pairs (recently emerged from 
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mitosis), Do(t) the number of old cell pairs, and T(t) the amount of toxin in contact 
with them. Then we write 







+He-W\C-D n (t)-D 0 (t)) 


(3.16) 


d °f=-X 0 e^<) m + to M 


(3.17) 


.. v„T(l)D r (t) v 0 T(t)D 0 (t) 
dt U 1 +K n T(t) l+K 0 T(t) 


(3.18) 



Consider first (3.16). The first term on the right-hand side represents a death 
rate of new cell pairs that increases with toxin concentration and is zero with 
none. It represents the effect of toxin-induced necrosis upon the new cell pairs. 
The second term is an aging term: it represents the rate at which new cells enter 
the old-cell population. The third term represents the mitosis rate; mitosis 
resembles birth in that new cell pairs are the product. It has been assumed that 
mitosis rate always decreases with increased toxin concentration, but such need 
not be the case. 

Next, consider (3.17). The first term is a death rate for old cell pairs that 
increases with toxin concentration; it represents both necrosis and apoptosis. The 
second term is the rate of maturation of new cells into old. Doubtless the 
transition rate, <p, could depend upon toxicity in the organ, but this possible 
dependence is not modeled here. Biological insight and information can perhaps 
suggest an appropriate and interesting functional dependence for <p. 

Finally, (3.18) models the rate of toxin concentration increase. The first term is 
the rate at which the basic toxin reaches the vicinity of the cells. The second and 
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third terms both represent the rates at which toxin concentration is reduced by 
action of new cells (second term) and old cells (third term). Note that both cell 
typ s reduce toxin concentration roughly in proportion to that concentration, but 
their effects saturate a la Michaelis-Menton. 

The equations supplied are illustrative of what may happen for certain toxins 
and cell types; they are intended to provoke discussion and stimulate 
experimental and observational checking and revision. 

A steady-state solution to equations (3.16) - (3.18) for a constant toxin input 
t it) = t would satisfy 



0 = - (e^ T - l)D„ - <PD» + lie-* 1 (C - D„ - D 0 ) 
0 = -A 0 ^ oT D 0 + <pD n 

0 = T v n TD n V0TD 0 

1 + XjjT 1 + Kq T 

Simplification of (3.16a) - (3.17a) results in 



(3.16a) 

(3.17a) 

(3.18a) 



Dn = —e~ ri ° 1 D r 



(3.19) 



D n =li 0 e S T C 



X n {e Tlnl -lj+ ip + /le 1 + j-e 7,07 



-l 



(3.20) 



Finally, expressions (3.19) - (3.20) can be substituted into equation (3.18a) to 
obtain an equation for T. This latter equation can have 0, 1 or 2 solutions. If it has 
no solutions, then D n = Dq = 0 and the organ is dead. If the equation has two 
solutions, then the smaller of the two possible solutions T s is biologically 
plausible since the smaller solution is an increasing function of r, whereas the 
larger solution is a decreasing function. 
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If r = 0, then T = 0 and 







(3.21) 



(3.22) 



4. Numerical Illustrations. 

We have submitted the above hypothetical equations to several dosage 
regimes. 

Illustration 1. 

Suppose toxin input is 0 until time t = 5, after which a constant exposure t is 
given. Figure 1 displays the function 



where D„ and Dq are obtained from equation (3.19) and (3.20); the parameter 
values used are displayed on the figure. The (possible) solution(s) to equation 
(3.18a) can be found drawing a horizontal line at the value of T of interest. The 
limiting amount of toxin in the organ subject to a constant input of toxin r is the 
smaller of the two values of T corresponding to the intersection points of /with 
the horizontal line. Note that there is no solution for r > 3.53, the maximum value 
of the function. This means that toxin has killed all old and new cells; long before 
that occurs the organism has died. 

Figures 2-5 display the results of a numerical solution of equation (3.16)— 
(3.18) for values of r = 0, 1, 2, 3, 3.3, and 3.53. The other parameter values used 



flj) — + tJpTDo 

n ' 1 + K n T \+k 0 T 
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are displayed on the Figures. The initial values of the numbers of young and old 
cell pairs and the amount of toxin present are the limiting values when x = 0. 

Figure 2 shows the effects of various levels of toxin on the young cell pair 
population. At the beginning of exposure the young cell population decreases. 
However, if the level of toxin is not too large, after this initial decrease the young 
cell population increases to a limiting value higher than that occurring when 
x = 0. If the level of exposure is too high, the young cell population decreases to 
zero. 

Figure 3 shows the effects of various levels of toxin on the old cell pair 
population. If the level of exposure is too high, the old cell population decreases 
to 0. For low or moderate exposure levels the old cell population initially 
decreases but then recovers somewhat to a limiting value which is below the 
limiting value if there were no toxin. 

Figure 4 displays the total number of totally differentiated cell pairs. As the 
level of exposure increases the total number of cell pairs decreases. The total 
number of totally differentiated cell pairs initially decreases, then recovers 
somewhat before decreasing a little again to a limiting value. 

Figure 5 displays the level of toxin in the organ. 

It seems possible that our model represents an effect called hormesis wherein 
a smallish amount of deleterious agent induces a biological entity to over- 
compensate for a (low-level) insult, but eventually to succumb to a larger dose. If 
young cells are actually more efficient and productive than old then the 
proportional buildup of the new-cell population for low toxin levels has just such 
an effect in our example. 
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5. Stochastic Versions of Models 

The purpose of this section is to describe a fully stochastic version of the cell- 
toxin interaction. Note that once stochastic elements are permitted to enter a 
great many optional behaviors can be modeled. It will be seen that expressions 
(3.12), (3.16), (3.17), and (3.18) are not exact replicas of the correct mean or 
expected -value equations. 

Model S.l. Markovian Model of Cell-Toxin Interaction. 

Suppose D(t ) represents the number of totally differentiated cell pairs in an 
organ with a maximum of C pairs, and T(t) is the toxin concentration at time t. 
Both D(0 and T(f) are now considered to vary randomly. 

Let (D(f), T(t), t > 0) be a bivariate birth and death Markov process in 
continuous time and with state space R x R, i.e., all pairs of integers such that 
0 < D(f) < C, 0 < T(t). Then put for the transition probabilities from / -* t + A 

P{(D(I + A),T(I + A)) = (D(t) + 1,T(/))[D(/),T(/)J = p (T(l))(C - D(f))A +o(A) 

P{(D(f + A ),T(I + A)) = (D(l) - l,r(f))|D(f),r(f)} = A(T( f))D(f)A + o(A) 

P{(D(t + A), T(l + A)) = (D(l), T(l)+ 1)|D(/),T(()} = r(l)A + o(A) 

P{(D« + A),T(t + A)) = (D(f),r(0-l)|D(/),r(0} = P^|jA + o(A) 

P{(D(f + A),T(/ + A)) = (D(1),T(/))|D(|),T(0} = 1 - p(D(t),T(t))A + o(A), (5 ]} 

where p(D(t),T(t)) =/l(T(l))(C-D(t)) + A(T(l))D(l) +r(l)+ vD(t)T(t) / (1 + kT(|)); 
1-pA is simply the probability of no change in the time period (f, t + A). 

Examination of (5.1) reveals that, under very special conditions, expected- 
value and higher-moment equations can be written down, and can be solved 
numerically. Non-linearities in general impede such a step, but approximation 
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and perturbation methodologies can be invoked to produce useful 
approximations. 

The Mean Equations 
Invoke (5.1) to obtain 



E[D(( + A)|D(/),T(0] = (D(f) + lMr(t)XC - D <0)A 

+(D(0-l)A(T(i))D(0A 

4 D(f)[l -/l(T(f)XC - D(0)A - ;.(T(/))D(/)A] + o(A) 

= p(r(0 )(C - D(0)A - A (T(/))D(1)A + D(f ) + o(A). (5 2) 



Hence, upon removing the condition and passing to the limit 



J t E[D(0] = £[/i(r(f ))]C - £[D(/)/i (r(/))] - £[;. (t(;))d(i)] 



Next, 

£[T(I + A)]D(l),T(l)] = (T(i)+ 1)<()A + (T(0 - 1) V A 



+r(0 



1 - 



v ' 1+ vT(f) 



y. 



Removal of the condition leads to 



(5.3) 



|'Elr(0]=r(()-^E 



mm 

1 +KT(t)_' 



(5.4) 



These resemble the previously-presented deterministic model (3.12), but only 
roughly, since expectations of non-linear functions such as appear above cannot 
be directly evaluated. Special methods can be used that give satisfactory 
approximations. 
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Clearly the equations (5.3) and (5.4) cannot be solved without the introduction of 
approximate equations for the non-linear expectations on the r h s. One approach 
for doing this is presented in Appendix A. A second approach is presented in the 
next section. 

Asymptotics for C Large. 

Since the number of cells in an organ, C, is large, being of order 10 7 , it is 
natural to consider explicit asymptotic approximations based on that fact. In 
what follows we proceed formally. For precedent here see McNeil and Schach 
(1973) but also more recent work. 

Define the joint moment generating function (assumed to exist), to be 



(d) r 

v{d d ,e t ;t) = E 



,e d D(t)+6J(t) 



(5.5) 



Use the Markov property to derive forward equations in transform space: 



e 6 d D(t+A)+6,T(t+A)lp^ ,T(f)j = 



e 0 d( D (O + l) + 0 i’ r (O^(j(^))(c _ D(f))A + ^ +e ‘ T ^?.(T{t))D[t ) A 



|c e,D(tM,(r(Q+i) T(m , Whe^no-i) vD(t)T(t)A 
W l+\T(f) 



+e Wt)+9 t T(t)^ _ )) A ] + 0(A), 

where 

p(D(t),T(t)) = n(T(t)){C - D(l)) + A (T(O)D(i) + f (f ) + 
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= '¥(e d ,e i ;t + A) 



Thus 



,e d D(t+A)+9 t T(t+ A) 



= £ 
+E 



e^D(0+^+W) /z (j(^)( C _D(f))A] 



,e d D(t)+e t T(t)-e t uP(Qr(Q A 
1 + jcT(f) 



■^D(o+e,T(f) + 0 (T(f ) A J + 
^D(0 + e ( T(0]_ E [ e ^D(/) + 6 i r(0 p(D(f)/T(f) )] 



Consequently, after re-arrangement and division by A, and letting A — » 0, 



= ( t <fc - ij^WW^c- D((»] 
-l)E[e e ‘ ,D! ' )+e ' T( 'U(T(())D(0] 

+(*»■ 



+ 



-i)e 



e rf p«)+e,r«) p (0 r (0 
l+»cT(f)_' 



(5.7) 



These equations are highly non-linear and intractable; consequently we introduce 
an appropriate asymptotic analysis that assumes C, the number of cell pairs, to be 
large. 



Scaling 

Define 



x( ,) = P(0-Ctt(0 y(i) - r(i)|M 



(5.8) 



and the joint mgf of the scaled "noise" variables X(f) and Y(t): 
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<p(8dA; 0 



(d) 



E 



e e d x(t)+e t Y(t) 



= E \ e Odm/'lc+mt)/Jc 



e -(e d Jca(t)+e t Jcp(t)) 



= v(e d / Vc, e t / 



(5.9) 



Consequently, 

x v(e d /Jc.ej Vc,/) = (p{e d ,d t -,t)e^ edai < l)+6 ^\ (5.10) 

It is necessary to find equations for a(t ) and P(t), and the joint mgf (p(6 d , 6f, t). 
Proceed as follows: 

use (5.7); re-define the transform variables as below 6 d -» 6 d / VC , 6 t -> 6 t / VC ; 
from (5.10) we get 




+ <pjc(e d a-(i)+ e t nt))‘^ (6M,) * e,m 







/Vc)(Vcx(0+ca(o) £ ,(^/Vc)(Vcr(0+c^(0) 



x{(> /Vf - i)(/i(c/)(0 + Vcr(i))(c - Co(0 - Vcx(f ))) 
_lV;.(C^(/) + ,,/cy(f))(Co(f)-VcX(/))) 




,- 9 , /Vc _ 1 \(c«{<) + ■Jcx(t))(cpQ)+ Vcv(0) ' ' 
1 + v(c/J(f)+Vcy( 0 ) 



(5.11) 



Next scale the transition rates 

|i(Q3(/) + VCY(O) = ^*(/?(f) + Y(f) / VC) (5.12a) 

A(Q3(f ) + VCY(O) = (£(/) + Y(f ) / VC) (5.12b) 
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(5.12c) 

(5.12d) 



v(t)=v*(t)/C r(t) = c/(t) 

K(cp(t) + VCY(O) = K(p{t ) + Y{t) / Vc) 

in order to obtain non-degenerate results involving all aspects of the process as C 
becomes large. Starred functions or parameters are 0(1). 

Expand state-dependent parameters by Taylor series in powers of 1 / VC after 
cancelling e 'JC{0 l ja{t)+9 t p{t)) f rom both sides; we put Uo{P(t)) ^ or the ^st term of 
the expansion of jz*(j3(f) + Y(f)/Vc); notation for X and k* is analogous. We 

obtain 



4£ + Wc(^a'(f)+W)) = 

e J x(/) c e,y(o x |^ + || + . ^;( W) ) + ^;( W) ).^..^ c _ Ca(() _Vcx(0) 

"Jc + 2c + " + mW) ■ + ^M) 



3t 

E 



e, e} A 

—7= + — - — h 

VC 2C 



Cr(f) 



0, 0f 
VC 2C 



u*(0/C 



(Ca(f) + VcX(f))(Cjg(f) + VCY(f)) 



1 + K o{P( t )) + *i(0(O) y (O/ 

Re-write the last line above in expanded form to get 

+ ]-’(0(«(0 - xco/ vc + . ,.)(pco- no / vc + .. - 

f 1 ” (*i (0(0) / ( 1+ K 'o(0(O)))^r + - • -1 



(5.13) 



(5.14) 
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&+-Jc(e d a :'(!)+ WO)? = 



(f 

2C 



c^of'K 1 «(0)« ,+ /‘i( , X 1 “WjcM' Uo ^Tjcdei 



1 d(D 



+. 



1l + ^} 

VC 2C 



c 



Ao( ( )a«)«- + A;(,)a(0^|| + A' o (,)-^^ + ..' 



e, ef ' 
-4= + — -+... 
Vc 2C 






0, 0? ' 

— p - + — +... 

VC 2C 



V (t) » w 0 + 

i+*o (0 



f(V(owl 


1 dtp ( 


l 1 + r o('h 


/ 



M 



+t ,- (() _ellL ’ is.. u - (l) £mihlM ’ i£. + .. 

i + K 0 (t)Jcde, (i + Ko(t)) ^ 3B < 



J) 



We express the dependence of <p on C as follows: 

(p(e d ,e t ;t,C)= (po{O d ,6 t ;t)+ (p ] (e d ,e t ;t)-j^+... 



\ 

y 



(5.15) 

(5.16) 
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Now substitute into (5.15) to obtain 



w((Xj- «(»))£ p/(^J 



;=o 

/ 

v 






VC 2C 



• M l v d( pj 1 



/ .i4. + fi + V 

VC 2C c 



lo(()o(0y — ■ — - + /. i (;)a(7)— ™ V — — 

0 ,to '(Vc)' ,w u Vc £,<*>, (vc) 






e t ef ' 

— v — ~ H K . . 

v VC 2C , 



Ct'Wlf/rprj- 

;=o (VC)' 



e, e, 2 ^ 

‘Vc + 2c + - 



c 



l+^o(0;VO '(VC)' 



+ 



»‘(0W) i y^i i 

i+^(0 -fc fc, se t (-Jc)’ 



+ 



^(Qg(0 
1+ *o(0 



1+ K Q (t) 



1 ^ <fy; 1 



M 



3c £„ 



“o ae < (Vc)' 



-+... 



where the omitted terms ("...") are 0(1/0, and Aq(0 = Aq(/?( 0) etc. 



(Vc)' 



(5.17) 
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Identify terms of order VC : 

lhs: (e d a'(t)+e,i3'(i) Wo 



h s: e d ^o(p{t))( 1 - or(0))- a{t)lo(P{t))](po + e t 



T *(0~ 



V(t)a(t) 0 (t) 

1+ *5(0(0) 



(5.18a) 
(Pq; (5.18b) 



since these must cancel, the functions exit) and (Kt ) satisfy the differential 
equations 



«'(0 = “ a ( 0 ) ~ A 0 W 0) a (0 






«W(0 

1 + K o[P(t)) 



(5.19) 



the latter must in general be solved numerically, but steady-state information can 
be derived more easily. The system of equations (5.19) is a scaled version of 
equation (3.12). 

Next, identify terms of order 1 (coefficient of (l / Vc) . The expression 
obtained is of the form 



(O+^-MO 



Si 



<P 0 






(5.20) 
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Here the above coefficients are 



Ai(0 = (^o(0(l-a(0) + ^o(0a(0) 



MO = *{*) + 



v{t)a{t)P{t) 



1+ K Q (t) 

Bdd(t) = -/*o(0-*o(0 

B ft) - v ‘(0i8(0 

l+K'o(0 

B/d(0 = «(0)-*i(0«(0 

B m = ^(f)g(o r p{t)K]{t) 

1 + K o(i) _ 1 + K o(0 . 



(5.21a) 

(5.21b) 

(5.21c) 

(5.21d) 

(5.21e) 

(5.21f) 



The resulting partial differential equation for <po is recognizable as characterizing 
the joint mgf of an Omstein-Uhlenbeck (Gaussian) stochastic process. 

Moment Equations. 

Differentiation of the pde with respect to 6d and 6t allows recovery of 
differential equations for the covariance function of (X(f), Y(0), the scaled 
stochastic terms describing deviations from the mean (to order C) term 



E[D(0] = Ca{i), E[T(/)] = Cp{t) (5.22) 

Differentiation of <Po(&d'&t''0 0^ = 6 t =0 shows first that if X(0) = Y(0), then 

E[X(t)l = E[Y(t)] = 0 to the order suggested. Second derivatives at 6 = 0 then 
deliver these differential equations for the variances and covariances. 
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Var[X(t)] — E[x 2 (f)] = m dd (t),Var[Y(t)] ^ E[Y 2 (f)] = m M (f) (5.23) 
Coy[X(0, Y(0] = E[X(0* Y(0] = m^f), (5.24) 

^ ) = A i (0 + 2B <fJ ( f ) m dd ( f ) + 2B td m dt (0 (5-25) 

= fc(0 +B //(0Kf(0+%(0^(0 +B ^i<(0 (5.26) 

^ (0 = A (0 ■ + 2B tt (t)m u [t) + 2 B dt m dt ( / ). (5.27) 



These equations can be solved numerically. Owing to the appearance of the 
Ornstein-Uhlenbeck, Feller (1966), we can argue that for large C, D(0 is 
approximately normal with mean Ca(t) and variance C m dd {t)-, T(t) is 

asymptotically normal with mean C Pit) and variance C mttit ); the two quantities 
are correlated with correlation coefficient p d {{t) = / ^nt dd (t)m tt (t) • 

Steady-State Solutions: Dose-Response for Steady Exposure 

Suppose t(t) = t , a constant, and that exposure has proceeded for some time 
so that a steady-state condition has been reached; this is modeled by letting 
o!it) = p’(t) = 0 in (5.19). We have, suppressing t and invoking the tentative 
working parametric forms (3.13), 






* ^(1 -a) = AqC 77 ^ a 



» * 
r = v 



aP 



(5.28a) 

(5.28b) 



l+K 0 (p) 

Notice that (5.28a) can be solved for p, the steady-state toxin concentration, in 
terms of a, the steady-state fraction of totally differentiated cells: 



P = TT^-rln 

£ +77 



Mol 


' 1 - a Yj 


Uo' 


v a )j 



(5.29) 



This can now be inserted into (5.28b) to obtain an implicit dose-response 
relationship: 
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*_ * a 1 

1+kJ(/J)^* + 7j’) Uov « J 



(5.30) 



the above, (5.30), can be solved for a[r*] so as to depict the fraction of totally 
differentiated cells' dependence on input toxin; the latter must here be a steady 
flow. 

Slope of Dose-Response Curve for Small Dose. 

Of interest to risk analysis is the behavior of the dose response curve for small 
values of (toxic) dose. We approximate this by finding an expression for 



da 



il 



d x 



_d_ 

t‘=0 dx 



rfl-CNTl) = , the rate of increase of the fraction of non- 

l l J / T * =0 



functioning totally differentiated cells when toxin input is very small. 

First notice that if x - 0 then Hq(\ - «) / AgCt = 1 . So differentiation of (5.30) at 
t* = 0 provides 



da 


* ' 
X 


f 


dx 


h’i 

o 



V v + 



(5.31) 



Another differentiation will give information concerning the curvatures of the 
response, or propensity to exhibit "hockey-stick" or knee-like or threshold 
behavior at small dose levels. A knowledge of such behavior is of interest to 
toxicologists and regulators who are concerned with risks in the workplace. It 
will be kept in mind that the behavior elucidated depends to an unknown degree 
upon the suitability of the models. 
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Model S.2. Markovian Model of Cell- Age-Specific Cell-Toxin Interaction. 

Suppose D n (t ) represents the number of young (or new) cell pairs (recently 
emerged from mitosis), D 0 (t) the number of old cell pairs, and T(f) is the toxin 
concentration at time t. All the variables D n (t), D 0 (t ) and T(f) are now considered 
to vary randomly. 

Let (D o (0, D n (t), T(t ); t > 0) be a trivariate birth and death process in 
continuous time and with state space R x R x R, i.e., all triples of integers such 
that 0 < D 0 (t ) + D n (t ) < C, 0 < T(t). Then, put for the transition probabilities from t 
t + A. 

P{(D 0 (t + A),D„(1 + A ),T(I + A)) = (D„(l) + 1,D„(/)- l,T(f))|D„(l),D„(/),T(f)] 

= 0 D„(l)A + o(A) 

P{(D 0 (l + A),D„(1 + A),T(1 + A)) = (D o (l)-l,D„(l),'T(l))|D ( ,(l),D„(0,T(i)} 

= X o (T(0)D o (0^ + o(A) 

P{(D o (f + A),D n (( + A),T(M-A)) = (D [ ,(0,D„(() + l,T(l))|D o (l),D„(l),r(0} 

))(c - D„ (f ) - D 0 (/))A + o( a) 

P{(D 0 (t + A),D„(( + A),T(f + A)) = (D„(0. D„(/) - 1,T(«))|D 0 (# ), D„ (f), T(f)} 



= A, ; (T(/))D„(1)A + 0(A) 

p{(D 0 (l + A),D„(/ + A),T(l + A)) = (D 0 (l),D„(l),T(() + l)!D 1 ,(l),D r ,(f).T(l)} 



= r(f)A +o(A) 

P((D 0 (I + A),D„(| + A),T(1 + A)) = (D o (l),D„(l),T(0 - l)|D„(l),D„(/),T(l)} 



U " T(I) D n (t) + T ^^ D 0 (I) 



A + o(A) 



l + K' n T(0“ nv ' / ■ l + v o T(0 
P{(D 0 (f + A),D„(l + A),T(/ + A)) = (D 0 (t),D„(l),T(f))jD < ,(f),D„(f),T(/)] 
= l-p(D 0 (l),D„(f),T(l))A + o(A) 
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where 



p(D o (l),D„((),T(f))=(>D n (l)+A o (T(0)D o (0+p(r(0)(C-D„(l)-D o (l)) 

+A„(r(/))D„ (()+<() 

+ rTm Dn(l)+ TT7Tm DM 

l+K n T{t) l + K 0 T(t) 

Asymptotics for C Large. 

In this subsection we invoke the techniques applied for the Markovian model 
of Cell-Toxin interaction to obtain a limiting trivariate normal distribution for the 
number of young and old cell pairs and amount of toxin in the organ as the 
number of cells in the organ C — > 

Define the joint moment generating function, (assumed to exist), to be 

{d) 

y{9 o ,d n> e t ;t) = E[exp{e o D o {t) + e n D n {t)+0 t T(t)}]. (5.32) 

Use the Markov property to derive the forward equations in transform space: 

E[exp{e 0 D 0 (f + A) + 0„D„(! + A) + 9,T(I + A)}|D 0 (<),D„(l),T(t)] 
=exp{eo(Do(i)+i)+e„(D„(f)-i)+e,r(0}«D„(i)A 
+exp{e <J (D 0 (f ) - 1) + 0„D„(f) + e,T(f)}^(T(0)D o (f)A 
+exp {e„D 0 (i) + e„(D„(f ) + 1) + 9,T(f)}p(T(/))(C - D„ (f) - D„(f )) A 
+exp{e 0 D 0 (t) + e„(D„(t)-i) + e t T(t)}^(T(t))D K (t)A 
+exp{e 0 D 0 (i)+e„D„(t)+e l (T(t)+i)]t(i)& 



+exp{e 0 D 0 (t)+e„D n (t)+ e,(T(i) - 1 )} 



v„mp„(t) v 0 T(t)D„(tj 

’L i+^no i+».' 0 t(o, 

+exp{e o D £ ,(l) + e„D„(f) + e f (T(l))][l-p(D IJ ((),D„((),r(/))A]+o(A). 



(3.33) 
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After rearrangement and division by A, letting A 0, and putting Z(t) = 
exp [6 0 D 0 Q) + 6 n D n (t ) + 0/7Yt)} 

dt 

= (e 6 ° - i)e[Z(/ )4>D n (/)] 

+(e- e » -l)E[Z(t)A o (7-(())D o (0] 

+(e e " - l)E[z(t)n(T(l))(C -D 0 (l)-D „ (/))] 
+(r e “-l)E[z(/)A„(T(/))D„(/)] 

+(* 9 ‘ -l)r(/)E[Z(i)] 



- i)e 



Z(0 



D,T(t)P, ; (Q , VJWM 



+ * 



l + v r T(/) l + v 0 T(/) 



Define as before 



(5.34) 



x o (0 = x„(/) = y(f) = 

VC VC VC 

and the joint mgf of the scaled "noise" variables X 0 (t), X n {t), Y(t): 

(d) 

v(e 0 ,e,,e i: t) = E[ex-p{e„x () (/)+ e„x„(/)+ e,r(0}] 



(5.35) 



= 'F 



e n 



,-^,-^lexp {-4c(e 0 a 0 {t) + e n a n (t)+e,p(t))}. 



.Vc'Vc'Vc J 



(5.36) 
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Proceeding as in the derivation of (5.11) 



~ exp{ VC (d 0 a 0 (t ) + 6 n a n (t) + 6 t p{t ))} 

n^(0X(O + ^(O + ^/? , (O)exp{Vc(^a o (/) + 0 n a w (O + 0<)3(f))} 



= £ 



exp{( 0 o / Vc)(VcX 0 (f) + Ca o (0) + ( 0 „ / Vc)(VcX n (f) + Ca H (t)) 



h(^/Vc)(Vcy(f)+c/?(f))} 






e fe-«,)/Vc_ 1 'j^ Can+V c X(i(() ) 



•o 0 /'lc _ i^„(c/3 (/) + Vcy(/))(ca 0 (/) + Vcx o (0) 

Y e " 1 ^ - {y(m) + -JCY(I))(C(1 - aM - «„(/)) - Vcx 0 (() - Vcx„(0) 
^ - i)*„(cW ) + VC y(f ))(ca„ ( 0 + Vex, (0) 

+fc 9 ' / '^-llr(l) 



r ^-e,/vc_,j 



«»(c/?(t) + >/cy(l)XCa n (<) + VcX„(()) ' 
1 + K 'n(Cfl(0 + '' l C Y ( 1 )) 

v c (Cftt) + JCY(l)'KCaM + &X < ,(t)) ' 
l + «r 0 (C/J(f) + Vcy(f)) 



(5.37) 
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Next scale the transition rates as before 



m(^P( ! ) + VCV(O) - MoWO) + Ml 


(5.38a) 


K( C M + Vcy(0) = a;, o (j3(0) + 4i(/J(i))^ 


(5.38b) 


A„(c J Q(0 + VCY(/)) = a;,o(/3(0) + v n/1 (/3(f))^ 


(5.38c) 


v = v /C 


r(0 = Cf*(f) 


(5.38d) 


u 

II 

V 


ii 

n 


(5.38e) 



Substituting into equation (5.37) and identifying terms results in the following 
equations. 

The terms of order Vc imply that a 0 (t), a n {t), and Pit) satisfy the differential 
equations 

«o(0 = ^ a n(0“^ooWO) a o(0 (5.39a) 

«n(0 = ^o(^(0)( 1 - « o {t)-a n (t))-(pa n {t)-?* n0 (p{t))a n (t) (5.39b) 

PV) = r(l)- - vj)(t)a n (t)- l— v'MM)- (5.39c) 

l + K n P{t) l+K 0 P{t) 

These equations are scaled versions of equations (3.16) - (3.18). 

The terms of order (Vc)° result in the following equations. 
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where 



*PQ = 

dt 



A>(0 + - ^ L Ai(0 + 2^(1) ^o^n < t >a n 



<Po 



) 



+e, 



+o. 



+ 0/ 



B_-^ + B„ „4^ + B, 






^_t x D 

o,< 



<90 



< J 



B 



f^O+B ^a + B 

n>0 ddn n,n do, 



d(p Q 

» ^ J 



Pt.i 






<?0n 



<,n 



dd r 



J t,t 



90 , 



A 0 { 0 — 0 < ^n(O'*’ ^o,o(0(O)®o(O 

A. (0 = ^ (j8(0)[l' - « 0 (O' - a„ (0] + A„ #0 ( 0(OK (O' + <t> a n (0 



A(0 = 






1 + K n P(t) 



0(O«n(O + 



t»n 



1 + *o0(O 



P{t)a 0 {t)+T 



\o=<M t )) 

Bo,n = 0 

B 0 ,t=-**o, l(0(O)«o(O 
Bn,o = So{m) 

B„ / n=-/ioWO)-4oWO)-0 

B n,f = /<l WOK 1 " «o(0“ «n(0)~ «n(0^n,l(^(0) 
_ -Up/?(Q _ -Uk0(O 

Of O ““ f tl "" r 

1 + ^(0 1 +*„«<) 



-Vn<*„( 0 


! *i»0(O 


v 0 a 0 (t) 


| * 00(0 


i+ *y3(f) 


\+K n p{t)_ 


1 + >c o 0(f) 


l+K" 0 /?(f)_ 



(5.40) 

(5.41a) 

(5.41b) 

(5.41c) 

(5.41d) 

(5.41e) 

(5.410 

(5.41g) 

(5.41h) 

(5.41i) 

(5.41 j) 
(5.41k) 



32 



Finally, calculations similar to those leading to (5.23) - (5.27) result in 



| E [ X »('>K 


+2B o . o e[x 2 (0] 

+2B 0 ,„E[X„(f)X n (/)] 






+2B o , J E[X„(t)T(0] 


(5.42a) 


^e[x 2 (/)] = .4„ 


42B„, o E(X o (/)X„(0] 

+2B n< „E[x 2 (l)] 






+2B„ -1 E[X„(f)T(f)] 


(5.42b) 




+2B ( , o E[X o (#)7-(0] 

+ 2B,,„E[X„(()T(0] 






+ 2B U E[T 2 (()] 


(5.42c) 



^E[X o (f)T(0] = [B„,„ + B u }E[X 0 (t)T(t)} 

+ B o ,„E[X„(l)T(l)] + B lj „E[X <) (l)X„(l)] + B„, 1 E[r 2 (0] 



+b,, 0 e[x 2 (i)] 



(5.42d) 



^■E[X„(0T(f)] = B n . o E[X o (f)T(f)] + (B„,„ + B t , ; )E[X„(0T(0] 



+B t , £ ,E[X o (f)X n (0]+B I ,„E[x2(()] + B„, f E[T 2 (0 



(5.42e) 
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f E[X o (0X„(0] = (B„ + B„,„ )E[X 0 (<)X„ (<)] 
+B„,E[X 0 (l)r(l)] 

+B„, ( E[X„(/)T(0] 

+b„,„e[x 2 (i)] 

+b 0 ,„e[x„ 2 (i)J 

—<pcc n 

The parameters in equations (5.19), (5.25) - (5.27) are as follows: 
For the case }i{T) = fJe'F, 

S(P (0 + Y{t) / VC) = H exp{-Z‘(p(t ) + Y(t) / VC)} 

where 

= £C 

Mo(P(t)) = 

Similarly, if A 0 (T) = ^Qe r, ° T , then 

^(0(0 + r (/) / VC ) = A„ expj rfe(P(0 + Y(l) / Vc)} 

= 4o(/3«)) + x;,,(p(l))^ 

where 

Vo ~ 7 lo^~ 

^.o(«0)=V* (,) 

^.iW0)=^v*w. 

Also 



(5.420 
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where 



A*„ (p{t) + Y(f) / VC ) = A n [exp{r]’ (jS(t) + Y(t) / Vc )} - 1] 

= Vo(£(0) + ^n,l(£(0)-^ 

*ln = T/ n C 

a;^,08(/)) = - 1“ 

4/l(/5(0) = ^n^ (O - 

K 0 = K' 0 C 

= *71^ 

t’(0 = t(0/C 
t/ = uC. 

A steady-state solution to equations (5.39a) - (5.39c) for a constant toxin input 
r(t) = r would satisfy 

0 = <pa n - X 0 e Ti °Pa 0 (5.43a) 

0 = fie~^ P{\- a 0 - a n )~ <pa n - -lja n (5.43b) 

0=t r- v n pa n v: v lP a 0 (5.43c) 

1 + K n P 1 + K 0 P 

where we are using functions of the form 

A„(7>A„(c’ , ” T -l) 

V(T) = ne-i T 

with 
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Simplification of (5.43a) and (5.43b) results in 



a o=J- e VoPa n 

A o 



(5.44a) 





-l)+<p + ne * P 


( A - n V 
l + ± e -Vj3\ 




v ) 


k A o )_ 



Finally expression (5.44a) and (5.44b) can be substituted into equation (5.43c) to 
obtain an equation for p. This latter equation can have 0, 1, or 2 solutions. If it has 
no solution, then a n = Oq- 0 and the organ is dead. If the equation has two 
solutions, then the smaller of the two possible solutions is biologically plausible 
since the smaller solution is an increasing function of r, whereas the larger 
solution is a decreasing function. 

Numerical Examples 

The following parameter values are used in the numerical examples below. 

Vo — 0.05, jJ. — 0.5, Tjo = 1, £ = 0.5, Kq = 1/ a o = 0.08, V n = 0.1, TJ n = 0.5, — 1, 

A„ = 0.05, d = 0.1, C = 10 7 . 



Figures 6-9 present results for the case in which xit) = 0 for 0 < f < 2 and then 
r{t) = 4 x 1CP for t > 2. The moments are started at their steady state values for no 
toxin input. Figure 6 presents the mean number of old cell pairs, C x o^it), and 

il/2 



the standard deviation of the number of old cell pairs. 



CxE 



(x 0 2 (<)]] • 



Note that 



the mean number of old cell pairs decreases to a new steady state value below 
the value for no toxin; the standard deviation increases, then decreases, then 
increases again to a new steady state value which is higher than the value for no 
toxin. Figure 7 presents the mean number of young/new cell pairs, C x a n (f), and 



the standard deviation of young/new cell pairs [c x E[x 2 (f)J] 1/2 . The mean 

number increases to a steady state value larger than that for no toxin input; the 
standard deviation initially decreases, then increases and finally decreases to a 
new steady state value below the value for no toxin. Figure 8 presents the mean 
number of totally differentiated cell pairs, C x (o^ff) + a n it)), and the standard 
deviation of the number of totally differentiated cell pairs 



[c x [e[x|(o] + 2E[x 0 (»)x n (o] + £[x„ 2 (f)]' 



-l1/2 



The mean number decreases to a new steady state value below the value for no 

toxin; the standard deviation increases to a new steady state value. Figure 9 

presents the mean amount of toxin, C x Pit), and the standard deviation 
r 2 11 1/2 

CxET (t) . Both values increase to a steady state value. 



Figures 10- 13 present results for a large input of toxin for a short time 
period: approximately a bolus input. The input of toxin is Tit) = 40 x 10 5 for 
2 < t < 3 and Tit) = 0 otherwise. Figure 10 presents the mean number of totally 
differentiated cell pairs and the standard deviation of totally differentiated cell 
pairs. The mean number initially decreases, then increases to the steady state 
value with no toxin. The standard deviation initially increases, then decreases to 
the steady state value with no toxin. Figure 11 displays the mean number of old 
cell pairs and the standard deviation of old cell pairs. The mean number initially 
decreases, then returns to the steady state value with no toxin. The standard 
deviation initially decreases, then increases, then returns to the steady state value 
with no toxin. Figure 12 displays the mean and standard deviation of the number 
of young /new cell pairs. The mean number initially decreases, then increases, 
after which it returns to the steady state value with no toxin. The standard 
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deviation exhibits similar behavior but on a different scale. Figure 13 presents the 
mean and standard deviation of the amount of toxin. 

Figures 14 - 17 present results for the same parameter values and same toxin 
input as Figures 10-13 except that (p = 0.01 instead of 0.1; that is, the mean time 
until a young/new cell pair becomes an old cell pair is 100 = 1/0.01 rather than 
10 = 1/0.1. This should increase the number of young/new cell pairs. Figure 16 
shows that the mean number of young /new cells is indeed larger. Figure 14 
presents the mean and standard deviation of the number of totally differentiated 
cell pairs. Note that the mean initially decreases, then increases to overshoot the 
steady state value with no toxin. The mean then decreases to the steady state 
value with no toxin. This behavior of overshooting the steady state value with no 
toxin has been observed in experimental studies; Portier (1993). 
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APPENDIX A 



Isham-VVhittle Theory 

In this Appendix we present another approach to approximating the non- 
linear expectations in equations (5.3) and (5.4). We follow the example of Isham 
(1991) and Whittle (1957) by (a) writing equations for second moments and joint 
transforms, and (b) "closing" the system of lower-moment equations by 
assuming that (D(t), T(t)) is bivariate Gaussian; we call this approach Gaussian 
closure (GC). The plausibility of this assumption follows from Model D.l above. 
Asymptotic perturbation methods in section 5 will show that a bivariate 
Gaussian process is of natural occurrence if the system of cells is large, as is true 
in practice. 

To proceed, we write second-moment equations: 



E 



D 2 



(< 



+ A)|D(f),T(0] = (D(0 + l) 2 /i(T(l)XC - D(0)A + (D(0- - if A (r(f))D(/)A 



+D 2 (()[l - (/i(T(/))(C - D(t))A - A(T(0)D(<) A)] + o(A). 



Simplification and limit-taking provides 



|e[d 2 (0] = 2E[D(f)/l(r(!))(C- D(0)] - 2E[d 2 (.-)2(T(0)] 
+E[#<(r(f))(C - D(l ))} + E[A(T(/))D(0] 



(A.l) 



e[t 2 (/ + A)|D(i),r(/)] = (r(f)+i) 2 t(/)a 



+(r(i)-i) 2 v p '^ r '') + t 2 (i) 

v ' 1 + KT(t) y ’ 



J <0a+ 3'Ma 

' v w 1 +kT(t) 



+ o(A) 
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from which 




(A. 2) 



Next 



E[D(( + A)T(f + A)|D(0, no] = (D(0 + l)(T(l)V(r(0XC - D(0)A 

+(D{1) - l)(T(l))A(r(/))D(l)A + D(f)(T(f) + l)r(f)A + D(,)(T(l)-l ) V ^^ | a 



+D(ono 

which leads to 



P(Qr(0 



1 -u(T(l))(C - D(1))A + A (T(l))D(l)A + 1(04 + u ] + ^ ( ' 



| E[D(f )T(0] = E[T(l)/i(r(<))(C - D(0)] - E[r(f)A (nO)D(O] 



+ T (t)E[D(t)\-vE 



D(0- 



mm 

1 + K-T(t) ' 



(A. 3) 



In order to go further, it is necessary to parameterize the rate parameters and 
the Michaelis-Menten (M.-M.) term. We set as before, for i, fi > 0, 



^(T) = /i 0 e^ r , A(T) = V rjr . (A. 4) 

If (D(f), T(f)) has a bivariate normal distribution, then the distribution is 
determined by the marginal means, the marginal variances and the covariance. 
Other moments can be expressed in terms of these qualities. For example 
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e[d 2 (0t(0] = e[d(i)] 2 e[7-(i)] 

+2E[D(l)]Cot'[D((),r(0] 

+ E[T(l)M D (0] (A .5, 

Equations (5.3) - (5.4), (A.1)-(A.3) can be evaluated using the following 
approximate numerical procedure. Put 

m D (t) = E[D(()], m T (t) = E[T(()], (/) = e[d 2 (()] m (!) = e[t 2 (I)], and 

ra DT (/) = E[D(f)T(()]. 

m D (t + A) = m D (l) + A{/i(m r (f ))C- (/l(m T (t)) + A(m T (f)))m D (l)} (5.3a) 



mj(t + D) = mj(t)+ A< r(t)~ v 



1 + Knijit) 



m 



dt(0 



(5.4a) 



m D 2i t + A ) = ™ D 2( t ) 

+A^2^(m T (t))C mD (t)-2^(m T (t)) + X{m T (t)fjm D 2 (f) 
+^i(m T (t))C + (A(m r (f))-/i(m r (f)))m D (f)} 



ni T 2 (f + A) = m T 2 ( t ) 



+ M 2 r()K(l) + m - 2 v t + J ■ | ) e[d(()T 2 (()] 
+ i + ™ T (0 mDT(0} 



(A. 2a) 



moj(t + A) = mpj(t) 



+A [n(mj (i ))Cnij(t) - (A (tn T (t)) + ^(wr(*))) w Dr(0 



(A. 3a) 



42 



The higher order moments E[D(0T^(f)] and E[D~{t)T{t )] are evaluated using the 
expressions for a bivariate normal distribution such as (A.5). 

Comparison of the Two Systems of Moment Equations. 

In this section we study the limit as C -» «*> of Equations (5.3) - (5.4) and 
(A.l) - (A. 3) using the scaling of the transition rates (5.12a) - (5.12d). The limits of 
the equations are as follows. 



^a(l) = ^(/J(l))(l-a(/))-Ao(«0)a(l) 



± p(l)=r -. v -jm!L 

dt HK> 1 +Kp(t) 



(5.3b) 

(5.4b) 



^£[x 2 (I)] = ^(/) + e[x 2 (I)]2B^ 

+2B,jE[X(f)y(0] 

-2a(0[pi' («())+ A', (/S(f))]E[X(0n0] 




(<)]=A(0 



+ 2B„E[y 2 (l)] 



+2Bj,E[X(l)y(l)] 



r i>V«(0ft( Q 

(l+K Pit))/' 




2 vfKt) 
1+^(0 



1- 



l+K p{t) 



E[X(t)Y(t)] 



(A. 2b) 
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